include <error.h>
include <gset.h>
include <imhdr.h>
include <mwset.h>
include <fset.h>
include <math.h>

# GETIMAGE files below.

include <pkg/gtools.h>
#include <smw.h>
#include <units.h>
include <ctype.h>
include <evexpr.h>      # For reading in data and evaluating the expression

include <clio.h>
include <imio.h>

task ldisplay

procedure ldisplay()

pointer immap(), imgl2r(), impl2r(), imps2r()

############################################################################
#
#  Variables for reading in from the .par file
#
############################################################################

define MAXLENS 1000

#
# Parametros
#
char input[SZ_FNAME]
char output[SZ_FNAME]
char ldispdir[SZ_FNAME]
char lconf[SZ_FNAME]
char lconfpath[SZ_FNAME]
real w1,w2

real crval1,cdelt1
int i1,i2

pointer fd
int  xl[MAXLENS]
int  yl[MAXLENS]
int  nl
int  xsize, ysize

pointer open (), fscan (), close ()
real clgetr()
real imgetr()
bool imaccf()

int stat

pointer ahivi()

pointer image_in
pointer image_out

int ndisp, nesp
pointer line
real linesum
pointer pixels

int i,j

pointer gp, gopen()

real sumrange()

char    command[SZ_LINE]        # Cursor command string
real    wx, wy                  # Cursor coordinates in WCS
int     clgcur()                                                           
int     wcs                     # Graphics WCS
int     key                     # Cursor key value

char clcmd[SZ_LINE]
char strfib[10]

begin

  call clgstr ("input", input, SZ_FNAME)
  call clgstr ("output", output, SZ_FNAME)
  call clgstr ("ldispdir", ldispdir, SZ_FNAME)
  call clgstr ("lconf",lconf, SZ_FNAME)

  w1 = clgetr("w1")
  w2 = clgetr("w2")

  call printf ("Parametros:\n %s\n %s\n %s\n %s\n %f,%f\n")
    call pargstr(input)
    call pargstr(output)
    call pargstr(ldispdir)
    call pargstr(lconf)
    call pargr(w1)
    call pargr(w2)

#  call malloc (imagecur, SZ_FNAME, TY_CHAR)
#  call strcpy ("imagecur", Memc[imagecur], SZ_FNAME)

  call strcpy(ldispdir,lconfpath,SZ_FNAME)
  call strcat("/",lconfpath,SZ_FNAME)
  call strcat(lconf,lconfpath,SZ_FNAME)
  fd = open (lconfpath, READ_ONLY, TEXT_FILE)
  nl = 0
  while ( fscan (fd) != EOF ) {
        call gargi (xl[nl])
        call gargi (yl[nl])
        nl = nl + 1
  }
  stat = close (fd)

  if (nl == 0) {
         call error (1, "Can't find any lens positions in the config file.\n")
  }

  xsize = ahivi(xl,nl) + 1
  ysize = ahivi(yl,nl) + 1

  call printf ("Found %d lens positions in the file '%s' \n")
                call pargi (nl)
                call pargstr (lconf)

  call printf ("IFU dimensions: %d x %d \n")
                call pargi (xsize)
                call pargi (ysize)

  image_in = immap(input,READ_ONLY,0)

  ndisp = IM_LEN(image_in, 1)
  nesp  = IM_LEN(image_in, 2)

  line = imgl2r(image_in, 16)
#  call gplotv(teste,9,1,9,"linha 16")
#  gp =  gopen ("stdgraph",NEW_FILE,STDGRAPH)
#  call gploto(gp,teste,9,1,9,"linha 16")
#  call gclose(gp)

  call printf ("Found %d spectra each with %d channels \n")
                call pargi (nesp)
                call pargi (ndisp)

  if (imaccf(image_in,"CRVAL1") && imaccf(image_in,"CDELT1")) {

    crval1 = imgetr(image_in,"CRVAL1")
    cdelt1 = imgetr(image_in,"CDELT1") 
    
    call printf("wavelenght calibration: w0 = %f ,dl = %f\n")
      call pargr(crval1)
      call pargr(cdelt1)
    call printf("Start wavelenght  = %f\n")
      call pargr(crval1)
    call printf("Ending wavelenght = %f\n")
      call pargr(crval1 + cdelt1 * (ndisp-1))

    call lambda_calib(w1,w2,crval1,cdelt1,ndisp,i1,i2)
 
  } else {
    call printf("Image not calibrated in wavelenght\n")
    i1 = 0
    i2 = ndisp - 1
  }

#  call printf ("Espectra no 16, channel 500: %f \n")
#                call pargr (Memr[line+500])

#  linesum = sumrange(line,0,ndisp-1)

#  call printf ("Espectra no 16, sum: %f \n")
#                call pargr (linesum)

#  image_out = immap(output,NEW_IMAGE,0)

  call display_image(i1,i2,image_in,output,xsize,ysize,xl,yl,nesp)
  call imunmap(image_in)

end

real procedure sumrange (array, lowval, highval)

pointer array
int     lowval, highval
real    sum
real    takeoff

int     i

begin
	if (lowval < 0)
	  return 0
        sum = 0.
        do i = lowval, highval
                sum = sum + Memr[array+i]
        return (sum)
end


procedure lambda_calib (w1,w2,crval1,cdelt1,ndisp,i1,i2)

real w1
real w2
real crval1
real cdelt1
int ndisp
int i1
int i2

begin

    if (w1 == INDEF)
      i1 = 0
    else {
      if (w1 < crval1) {
        call printf("Start wavelenght out of range\n")
	i1 = -1
        i2 = -1
	return
      }
      i1 = (w1-crval1) / cdelt1
    }

    if (w2 == INDEF)
      i2 = ndisp - 1
    else {
      if (w2 > crval1+cdelt1*ndisp) {
        call printf("Ending wavelenght out of range\n")
 	i1 = -1
	i2 = -1
	return
      }
      i2 = (w2-crval1) / cdelt1
    }
    if (w1 != INDEF && w2 != INDEF && w1 > w2) {
	call printf("Invalid input: w1 > w2\n")
 	i1 = -1
	i2 = -1
	return
    }

end

procedure display_image(i1,i2,image_in,output,xsize,ysize,xl,yl,nesp)

int i1,i2
pointer image_in
char output[SZ_FNAME]
int xsize,ysize
int xl[nesp]
int yl[nesp]
int nesp
	 
#real vin[IM_MAXDIM]
#real vout[IM_MAXDIM]

pointer image_out
pointer pixels
#pointer ptr1
int i
long one_l

pointer immap(), imgl2r(), impl2r(), imps2r(), impgsr()
real sumrange()
char cldisp[SZ_LINE]
bool imaccess()

begin

  if (imaccess(output)) {
    call imdelete(output)
  } 
  image_out = immap(output,NEW_IMAGE,0)
  IM_NDIM(image_out) = 2
  IM_LEN(image_out, 1) = xsize
  IM_LEN(image_out, 2) = ysize
  IM_PIXTYPE(image_out) = TY_REAL
#  one_l = 1
#  call amovkl(one_l,vin,IM_MAXDIM)
#  call amovkl(one_l,vout,IM_MAXDIM)
#  ptr1 = impgsr(image_out,vin,vout,2)
#  call imunmap(image_out)
#  image_out = immap(output,READ_WRITE,0)
  

  call amovkr (0.0, Memr[imps2r (image_out, 1, xsize, 1, ysize)], (xsize * ysize))

  pixels = imps2r (image_out, 1, xsize, 1, ysize)
  

#  Memr[pixels + 16+32*8] = 2.0

#  Memr[impl2r(image_out,1)+16] = 2.0


  do i = 0, nesp-1
     Memr[pixels+xl[i]+yl[i]*xsize] = sumrange(imgl2r(image_in, i+1),i1,i2)

  call imunmap(image_out)
  call strcpy("display ",cldisp,SZ_LINE)
  call strcat(output,cldisp,SZ_LINE)
  call strcat(" 1",cldisp,SZ_LINE)
  call printf ("disp: Executing the cl command:'%g'\n")
               call pargstr (cldisp)
  call clcmdw(cldisp)

end

